#ifndef CONVFFT
#define CONVFFT
#include <string>
#include <vector>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <algorithm>
#include "fft_gpu.cu"
#include "utils.cuh"
// Include Thrust
#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <thrust/transform.h>
#include<ctime>

/*
    a funtion to calculate the convolution of vector A and B 
    double *ina: vector A
    size_t N_A,: size of vector A
    double *inb: vector B
    size_t N_B : size of vector B
    double *out: vector out
    size_t N   : size of vector out
    size_t mode: 0 -- full, 1 -- same, 2 -- valid
*/
void convfft_1d(Real *ina, size_t N_A, Real *inb, size_t N_B, Real *out, size_t N, size_t mode)
{
    // calculate the length 
    // FFT_N = 2 ** (int(np.log2(YN)) + 1)
    // size_t L = pow( 2, static_cast<int>(log2( static_cast<double>(N) ) + 1 ));
    size_t L = N;
    // transform double to complex
    Comp* inaComp = (Comp*) malloc(L*sizeof(Comp));
    Comp* inbComp = (Comp*) malloc(L*sizeof(Comp));
    // init inComp as 0 vector
    memset(inaComp,0,L*sizeof(Comp));
    memset(inbComp,0,L*sizeof(Comp));
    transReal2Comp(ina,inaComp,N_A);
    transReal2Comp(inb,inbComp,N_B);
    // malloc output 
    Comp* outaComp = (Comp*) malloc(L*sizeof(Comp));
    Comp* outbComp = (Comp*) malloc(L*sizeof(Comp));

    Comp* incComp = (Comp*) malloc(L*sizeof(Comp));
    Comp* outcComp = (Comp*) malloc(L*sizeof(Comp));
    // give a and b a fft transform

    fft_cuda_1d(1,L,inaComp,outaComp); 
    fft_cuda_1d(1,L,inbComp,outbComp); 

    // multiply a and b to get c
    vectorCrossMultip(outaComp,outbComp,incComp,L);
    // give c a ifft transform
    fft_cuda_1d(0,L,incComp,outcComp);
    // fft_cuda_1d(0,L,outaComp,outcComp);
    // transform complex to double
    transComp2Real(outcComp,out,N);
}



/*
    a funtion to calculate the convolution of 2D vector A and B 
    
*/
void convfft_2d(Real* input_X, int in_n, int in_c, int in_d, int in_h, int in_w,
    Real* input_W, int k_n, int k_c, int k_d, int k_h, int k_w,
    int p_h, int p_w,
    int s_h, int s_w,
    int d_h, int d_w,
    Real* output_Y,int out_n, int out_c, int out_d, int out_h, int out_w)
{
 
    /** comments from cuDNN-API

      限制输入矩阵的通道数和kernel的输入通道数一致
      1. in_c must equals to k_c
      2. in_h + 2 * padding height must equal 256 or less
         in_w + 2 * padding width must equal 256 or less
      3. stride of height and width must equal 1
      4. k_h must greater then padding height
         k_w must greater then padding width

    */
    if (k_c != in_c 
      || in_h + 2 * p_h > 256 || in_w + 2 * p_w > 256 
      || s_h != 1 || s_w != 1 
      || k_h <= p_h || k_w <= p_w ){
      throw string("convfft_2d cannot perform convolution because the input size does not match the kernel size!");
    }
   
    // tmp size of output
    size_t tmp_out_h = in_h + k_h - 1;
    size_t tmp_out_w = in_w + k_w - 1;

    // allocate memory for the nchw output
    Comp* knkcOutOutComp = (Comp*) malloc( in_n * k_n * k_c * tmp_out_h * tmp_out_w * sizeof(Comp));
    // allocate memory for the nhw output ,because after accumulate to the first dimension, the 
    // k_c-th dimension will eliminate.
    Comp* knOutOutComp = (Comp*) malloc( in_n * k_n * tmp_out_h * tmp_out_w * sizeof(Comp));
    // malloc temp pointer
    // init temp in and kernel vector
    Comp*  inComp = (Comp*) malloc(tmp_out_h * tmp_out_w * sizeof(Comp));
    Comp*  keComp = (Comp*) malloc(tmp_out_h * tmp_out_w * sizeof(Comp));
    // malloc output 
    Comp* outinComp = (Comp*) malloc(tmp_out_h * tmp_out_w * sizeof(Comp));
    Comp* outkeComp = (Comp*) malloc(tmp_out_h * tmp_out_w * sizeof(Comp));
    // tmp in comp pointer
    Comp* inOutComp = (Comp*) malloc(tmp_out_h * tmp_out_w * sizeof(Comp));
    
    for (int nn = 0 ; nn < in_n ; nn++){
        int n_kn_kc_out_offset = nn * k_n * k_c * tmp_out_h * tmp_out_w;
        // the output channel of kernel filter
        for (int k_index = 0 ; k_index < k_n ; k_index++){ 
            int tmp_out_kn_offset = n_kn_kc_out_offset + k_index * k_c * tmp_out_h * tmp_out_w;
            int out_kn_offset = nn * k_n * tmp_out_h * tmp_out_w + k_index * tmp_out_h * tmp_out_w;
          // the input channel of kernel filter
            for (int c_index = 0 ; c_index < k_c ; c_index++){
                // the offset of each feature map
                int tmp_out_offset =  tmp_out_kn_offset + c_index * tmp_out_h * tmp_out_w;
                int input_offset = nn * in_c * in_h * in_w + c_index * in_h * in_w;
                int filter_offset = k_index * k_c * k_h * k_w + c_index * k_h * k_w;
                // set zero value for all elements in temp vector
               memset(inComp,0,tmp_out_h * tmp_out_w * sizeof(Comp));
               memset(keComp,0,tmp_out_h * tmp_out_w * sizeof(Comp));
                // set real value as complex  
                transReal2Comp(input_X + input_offset , inComp,tmp_out_h,tmp_out_w,in_h,in_w);
                transReal2Comp(input_W + filter_offset ,keComp,tmp_out_h,tmp_out_w,k_h,k_w);

                fft_cuda_2d(1, tmp_out_h, tmp_out_w, inComp, outinComp);
                fft_cuda_2d(1, tmp_out_h, tmp_out_w, keComp, outkeComp);

                // multiply a and b to get c
                vectorCrossMultip(outinComp, outkeComp, inOutComp, tmp_out_h * tmp_out_w);

                // give out vector a ifft transform
                fft_cuda_2d(0, tmp_out_h, tmp_out_w, inOutComp, knkcOutOutComp + tmp_out_offset );
        
            }

            matrixAddition(knkcOutOutComp + tmp_out_kn_offset , k_c, tmp_out_h, tmp_out_w);
            eliminateInputChannel(knkcOutOutComp + tmp_out_kn_offset, knOutOutComp + out_kn_offset, tmp_out_h, tmp_out_w);
        }


  }
    // 对输出结果进行裁边处理
    Comp* outOutTrimComp = (Comp*) malloc(out_n * out_c * out_h * out_w * sizeof(Comp));
    matrixTrimPadding(knOutOutComp, outOutTrimComp, out_n, out_c, tmp_out_h, tmp_out_w, out_h, out_w);
    
    // transform complex to real
    transComp2Real(outOutTrimComp, output_Y, out_n * out_c * out_h * out_w);
}
 /*
   make the compute processrun in a serial mode  
 */
void convfft_2d_v2(Real* input_X, int in_n, int in_c, int in_d, int in_h, int in_w,
    Real* input_W, int k_n, int k_c, int k_d, int k_h, int k_w,
    int p_h, int p_w,
    int s_h, int s_w,
    int d_h, int d_w,
    Real* output_Y,int out_n, int out_c, int out_d, int out_h, int out_w)
{
 
    /** comments from cuDNN-API

      限制输入矩阵的通道数和kernel的输入通道数一致
      1. in_c must equals to k_c
      2. in_h + 2 * padding height must equal 256 or less
         in_w + 2 * padding width must equal 256 or less
      3. stride of height and width must equal 1
      4. k_h must greater then padding height
         k_w must greater then padding width

    */
    if (k_c != in_c 
      || in_h + 2 * p_h > 256 || in_w + 2 * p_w > 256 
      || s_h != 1 || s_w != 1 
      || k_h <= p_h || k_w <= p_w ){
	    cout << k_c <<" " << in_c <<" " << in_h <<" " << p_h <<" " << in_w<<" "  << p_w <<" " << s_h <<" " << s_w <<" " << k_h <<" " << p_h <<" " << k_w <<" " << p_w <<" " << endl;
      cout << "convfft_2d cannot perform convolution because the input size does not match the kernel size!"<< endl;
      throw string("convfft_2d cannot perform convolution because the input size does not match the kernel size!");
    }
  
    // tmp size of output
    size_t tmp_out_h_ori = in_h + k_h - 1;
    size_t tmp_out_w_ori  = in_w + k_w - 1;

    // the size of h and w need to satisfy the power of 2
    size_t tmp_out_h = pow(2,((int)log2(tmp_out_h_ori - 1) + 1)); 
    size_t tmp_out_w = pow(2,((int)log2(tmp_out_w_ori - 1) + 1)); 

    // allocate memory for the nchw output
    Comp* knkcOutOutComp = (Comp*) malloc( in_n * k_n * k_c * tmp_out_h * tmp_out_w * sizeof(Comp));
    // allocate memory for the nhw output ,because after accumulate to the first dimension, the 
    // k_c-th dimension will eliminate.
    Comp* knOutOutComp = (Comp*) malloc( in_n * k_n * tmp_out_h * tmp_out_w * sizeof(Comp));
    // malloc temp pointer
    // init temp in and kernel vector
    Comp*  inComp = (Comp*) malloc(in_n * in_c * tmp_out_h * tmp_out_w * sizeof(Comp));
    Comp*  keComp = (Comp*) malloc(k_n * k_c * tmp_out_h * tmp_out_w * sizeof(Comp));
    // malloc output 
    Comp* outinComp = (Comp*) malloc(in_n * in_c * tmp_out_h * tmp_out_w * sizeof(Comp));
    Comp* outkeComp = (Comp*) malloc(k_n * k_c * tmp_out_h * tmp_out_w * sizeof(Comp));
    // tmp in comp pointer
    Comp* inOutComp = (Comp*) malloc(in_n * in_c * k_n * tmp_out_h * tmp_out_w * sizeof(Comp));
   


    // init the incomp and kecomp
    clock_t startTime,endTime;
    startTime = clock(); 
    transReal2CompWithMemset(input_X , inComp, in_n, in_c, tmp_out_h, tmp_out_w, in_h, in_w);
    transReal2CompWithMemset(input_W , keComp, k_n, k_c, tmp_out_h, tmp_out_w, k_h, k_w);
    endTime = clock();
    cout << "The transReal2CompWithMemset run time is: " <<(double) 1000 * (endTime - startTime) / CLOCKS_PER_SEC << "ms" << endl;

    std::cout << "tmp_out_h: "<< tmp_out_h << std::endl;  
    std::cout << "tmp_out_w: "<< tmp_out_w << std::endl; 

   cout << "keComp real part is " << endl;
   for (int nn = 0 ; nn < k_n ; nn++){
            for (int c_index = 0 ; c_index < k_c ; c_index++){
                // the offset of each feature map
               for (int i = 0 ; i < tmp_out_h ; i++){
		   for (int j = 0 ; j < tmp_out_w ; j++){
		       cout << (keComp + nn * k_c * tmp_out_h * tmp_out_w + c_index* tmp_out_h * tmp_out_w )[i * tmp_out_w + j].real()<<",";
		   }
		   cout << endl;
		}

		   cout << endl;
	    }
		   cout << endl;
    }
		   cout << endl;
    startTime = clock(); 
    cfft_plan plan_h = make_cfft_plan (tmp_out_h);
    cfft_plan plan_w = make_cfft_plan (tmp_out_w);
    endTime = clock();
    cout << "The make cfft plan run time is: " <<(double) 1000 * (endTime - startTime) / CLOCKS_PER_SEC << "ms" << endl;
    startTime = clock(); 
    // fft the incomp and kecomp
    for (int i = 0 ; i < in_n ; i++){
	for (int j = 0 ; j < in_c ; j++){
	    int offset = i * in_c * tmp_out_h * tmp_out_w + j * tmp_out_h * tmp_out_w ; 
            fft_cuda_2d_v2(plan_h, plan_w, 1, tmp_out_h, tmp_out_w, inComp + offset , outinComp + offset);
	}
    }

    for (int i = 0 ; i < k_n ; i++){
	for (int j = 0 ; j < k_c ; j++){
	    int offset = i * k_c * tmp_out_h * tmp_out_w + j * tmp_out_h * tmp_out_w  ;
            fft_cuda_2d_v2(plan_h, plan_w, 1, tmp_out_h, tmp_out_w, keComp + offset , outkeComp + offset);
	}
    }
    
    endTime = clock();
    cout << "The fft_cuda_2d run time is: " <<(double) 1000 * (endTime - startTime) / CLOCKS_PER_SEC << "ms" << endl;
    // vector cross multiply
    startTime = clock(); 
    for (int nn = 0 ; nn < in_n ; nn++){
        for (int k_index = 0 ; k_index < k_n ; k_index++){ 
            for (int c_index = 0 ; c_index < k_c ; c_index++){
                // the offset of each feature map
                int input_offset = nn * in_c * tmp_out_h * tmp_out_w + c_index * tmp_out_h * tmp_out_w;
                int filter_offset = k_index * k_c * tmp_out_h * tmp_out_w + c_index * tmp_out_h * tmp_out_w;

		cout << "input offset: " << input_offset << endl; 
		cout << "filter_offset : " << filter_offset << endl; 
                vectorCrossMultip(outinComp + input_offset, outkeComp + filter_offset, inOutComp + nn * k_n * k_c * tmp_out_h * tmp_out_w + k_index * k_c * tmp_out_h * tmp_out_w + c_index* tmp_out_h * tmp_out_w, tmp_out_h * tmp_out_w);
	    }
	}
    }
    endTime = clock();
    cout << "The vectorCrossMultip run time is: " <<(double) 1000 * (endTime - startTime) / CLOCKS_PER_SEC << "ms" << endl;
   cout << "outinComp real part is " << endl;
   for (int nn = 0 ; nn < in_n ; nn++){
            for (int c_index = 0 ; c_index < in_c ; c_index++){
                // the offset of each feature map
               for (int i = 0 ; i < tmp_out_h ; i++){
		   for (int j = 0 ; j < tmp_out_w ; j++){
		       cout << (outinComp + nn * in_c * tmp_out_h * tmp_out_w + c_index* tmp_out_h * tmp_out_w )[i * tmp_out_w + j].real()<<",";
		   }
		   cout << endl;
		}

		   cout << endl;
	    }
		   cout << endl;
    }
		   cout << endl;
   cout << "outinComp imag part is " << endl;
   for (int nn = 0 ; nn < in_n ; nn++){
            for (int c_index = 0 ; c_index < in_c ; c_index++){
                // the offset of each feature map
               for (int i = 0 ; i < tmp_out_h ; i++){
		   for (int j = 0 ; j < tmp_out_w ; j++){
		       cout << (outinComp + nn * in_c * tmp_out_h * tmp_out_w + c_index* tmp_out_h * tmp_out_w )[i * tmp_out_w + j].imag()<<",";
		   }
		   cout << endl;
		}

		   cout << endl;
	    }
		   cout << endl;
    }
		   cout << endl;
   cout << "outweightComp real part is " << endl;
   for (int nn = 0 ; nn < k_n ; nn++){
            for (int c_index = 0 ; c_index < k_c ; c_index++){
                // the offset of each feature map
               for (int i = 0 ; i < tmp_out_h ; i++){
		   for (int j = 0 ; j < tmp_out_w ; j++){
		       cout << (outkeComp + nn * k_c * tmp_out_h * tmp_out_w + c_index* tmp_out_h * tmp_out_w )[i * tmp_out_w + j].real()<<",";
		   }
		   cout << endl;
		}

		   cout << endl;
	    }
		   cout << endl;
    }
		   cout << endl;
   cout << "outweightComp imag part is " << endl;
   for (int nn = 0 ; nn < k_n ; nn++){
            for (int c_index = 0 ; c_index < k_c ; c_index++){
                // the offset of each feature map
               for (int i = 0 ; i < tmp_out_h ; i++){
		   for (int j = 0 ; j < tmp_out_w ; j++){
		       cout << (outkeComp + nn * k_c * tmp_out_h * tmp_out_w + c_index* tmp_out_h * tmp_out_w )[i * tmp_out_w + j].imag()<<",";
		   }
		   cout << endl;
		}

		   cout << endl;
	    }
		   cout << endl;
    }
		   cout << endl;

   cout << "inOutComp real part is " << endl;
   for (int nn = 0 ; nn < in_n ; nn++){
        for (int k_index = 0 ; k_index < k_n ; k_index++){ 
            for (int c_index = 0 ; c_index < k_c ; c_index++){
                // the offset of each feature map
               for (int i = 0 ; i < tmp_out_h ; i++){
		   for (int j = 0 ; j < tmp_out_w ; j++){
		       cout << (inOutComp+ nn * k_n * k_c * tmp_out_h * tmp_out_w + k_index * k_c * tmp_out_h * tmp_out_w + c_index* tmp_out_h * tmp_out_w )[i * tmp_out_w + j].real()<<",";
		   }
		   cout << endl;
		}

		   cout << endl;
	    }
		   cout << endl;
	}
		   cout << endl;
    }
		   cout << endl;

   cout << "inOutComp imag part is " << endl;
   for (int nn = 0 ; nn < in_n ; nn++){
        for (int k_index = 0 ; k_index < k_n ; k_index++){ 
            for (int c_index = 0 ; c_index < k_c ; c_index++){
                // the offset of each feature map
               for (int i = 0 ; i < tmp_out_h ; i++){
		   for (int j = 0 ; j < tmp_out_w ; j++){
		       cout << (inOutComp+ nn * k_n * k_c * tmp_out_h * tmp_out_w + k_index * k_c * tmp_out_h * tmp_out_w + c_index* tmp_out_h * tmp_out_w )[i * tmp_out_w + j].imag()<<",";
		   }
		   cout << endl;
		}

		   cout << endl;
	    }
		   cout << endl;
	}
		   cout << endl;
    }
		   cout << endl;

    startTime = clock(); 
    // fft 2d inverse
    for (int nn = 0 ; nn < in_n ; nn++){
        for (int k_index = 0 ; k_index < k_n ; k_index++){ 
            for (int c_index = 0 ; c_index < k_c ; c_index++){
                // the offset of each feature map
                fft_cuda_2d_v2(plan_h, plan_w, 0, tmp_out_h, tmp_out_w, 
				inOutComp+ nn * k_n * k_c * tmp_out_h * tmp_out_w + k_index * k_c * tmp_out_h * tmp_out_w + c_index* tmp_out_h * tmp_out_w, 
				knkcOutOutComp + nn * k_n * k_c * tmp_out_h * tmp_out_w + k_index * k_c * tmp_out_h * tmp_out_w + c_index* tmp_out_h * tmp_out_w );
	    }
	}
    }


    // destory the plan
    destroy_cfft_plan (plan_h);
    destroy_cfft_plan (plan_w);

    endTime = clock();
    cout << "The fft_cuda_2d inverse run time is: " <<(double) 1000 * (endTime - startTime) / CLOCKS_PER_SEC << "ms" << endl;
    startTime = clock(); 

   for (int nn = 0 ; nn < in_n ; nn++){
        int n_kn_kc_out_offset = nn * k_n * k_c * tmp_out_h * tmp_out_w;
        // the output channel of kernel filter
        for (int k_index = 0 ; k_index < k_n ; k_index++){ 
            int tmp_out_kn_offset = n_kn_kc_out_offset + k_index * k_c * tmp_out_h * tmp_out_w;
            int out_kn_offset = nn * k_n * tmp_out_h * tmp_out_w + k_index * tmp_out_h * tmp_out_w;
            matrixAddition(knkcOutOutComp + tmp_out_kn_offset , k_c, tmp_out_h, tmp_out_w);
            eliminateInputChannel(knkcOutOutComp + tmp_out_kn_offset, knOutOutComp + out_kn_offset, tmp_out_h, tmp_out_w);
        }
    }
    endTime = clock();
    cout << "The matrixAddition run time is: " <<(double) 1000 * (endTime - startTime) / CLOCKS_PER_SEC << "ms" << endl;
    startTime = clock(); 

  
    // 对输出结果进行裁边处理
    Comp* outOutTrimComp = (Comp*) malloc(out_n * out_c * out_h * out_w * sizeof(Comp));
    matrixTrimPadding(knOutOutComp, outOutTrimComp, out_n, out_c, tmp_out_h, tmp_out_w, out_h, out_w, tmp_out_h_ori, tmp_out_w_ori, 1);
    // transform complex to real
    transComp2Real(outOutTrimComp, output_Y, out_n * out_c * out_h * out_w);
    endTime = clock();
    cout << "The matrixTrimPadding and transreal run time is: " <<(double) 1000 * (endTime - startTime) / CLOCKS_PER_SEC << "ms" << endl;
    startTime = clock(); 

}








void convfft_3d(Real* input_X, int in_n, int in_c, int in_d, int in_h, int in_w,
  Real* input_W, int k_n, int k_c, int k_d, int k_h, int k_w,
  int p_d, int p_h, int p_w,
  int s_h, int s_w,
  int d_h, int d_w,
  Real* output_Y,int out_n, int out_c, int out_d, int out_h, int out_w)
{
    /**
      限制输入矩阵的通道数和kernel的输入通道数一致
    */
    if (k_c != in_c){
      throw string("convfft_3d cannot perform convolution because the input size does not match the kernel size!");
    }
   
    // tmp size of output
    size_t tmp_out_h = in_h + p_h * 2;
    
    // size_t tmp_out_d = in_d + k_d - 1;
    size_t tmp_out_d = in_d + p_d * 2;
    size_t p_out_w = in_w + p_w * 2;

    size_t tmp_out_w = p_out_w;
    if (p_out_w % 2 != 0){
      tmp_out_w = p_out_w + 1;
    }
    std::cout << "tmp_out_h: "<< tmp_out_h << std::endl;  
    std::cout << "tmp_out_w: "<< tmp_out_w << std::endl; 
    std::cout << "p_out_w: "<< p_out_w << std::endl; 
    std::cout << "in_w: "<< in_w << std::endl; 
    std::cout << "tmp_out_d: "<< tmp_out_d << std::endl; 

    if ((tmp_out_h - out_h) % 2 != 0 
        || (tmp_out_w - out_w) % 2 != 0
      || (tmp_out_d - out_d) % 2 != 0 ){ 
          // throw "output size is unsuitable, Please check the size of tensor!";
    }

    Comp* knkcOutOutComp = (Comp*) malloc( in_n * k_n * k_c * tmp_out_d * tmp_out_h * tmp_out_w * sizeof(Comp));
    Comp* knOutOutComp = (Comp*) malloc( in_n * k_n * tmp_out_d * tmp_out_h * tmp_out_w * sizeof(Comp));

    Comp* prePaddingComp = (Comp*) malloc(tmp_out_d * tmp_out_h * p_out_w * sizeof(Comp));
    Comp* inComp = (Comp*) malloc(tmp_out_d * tmp_out_h * tmp_out_w * sizeof(Comp));
    Comp* keComp = (Comp*) malloc(tmp_out_d * tmp_out_h * tmp_out_w * sizeof(Comp));

    Comp* outinComp = (Comp*) malloc(tmp_out_d * tmp_out_h * tmp_out_w * sizeof(Comp));
    Comp* outkeComp = (Comp*) malloc(tmp_out_d * tmp_out_h * tmp_out_w * sizeof(Comp));
    Comp* inOutComp = (Comp*) malloc(tmp_out_d * tmp_out_h * tmp_out_w * sizeof(Comp));
    // nckdhw
    for (int n = 0 ; n < out_n ; n++){// n 
      int n_kn_kc_out_offset = n * k_n * k_c * tmp_out_d * tmp_out_h * tmp_out_w;
      for (int k = 0 ; k < k_n ; k++ ){ // k_out

        int tmp_out_kn_offset = n_kn_kc_out_offset + k * k_c * tmp_out_d * tmp_out_h * tmp_out_w;
        int out_kn_offset = n * k_n * tmp_out_d *  tmp_out_h * tmp_out_w + k * tmp_out_d * tmp_out_h * tmp_out_w;

        for (int c = 0 ; c < k_c ; c++ ){ // k_in
          size_t tmp_out_offset =  tmp_out_kn_offset + c * tmp_out_d * tmp_out_h * tmp_out_w;
          size_t input_offset = n * in_c * in_d * in_h * in_w + c * in_d * in_h * in_w;
          size_t kernel_offset = k * k_c * k_d * k_h * k_w + c * k_d * k_h * k_w;



          memset(inComp,0,tmp_out_d * tmp_out_h * tmp_out_w * sizeof(Comp));
          memset(keComp,0,tmp_out_d * tmp_out_h * tmp_out_w * sizeof(Comp));
          memset(prePaddingComp,0,tmp_out_d * tmp_out_h * p_out_w * sizeof(Comp));
          
          inputPadding(input_X + input_offset, prePaddingComp, in_d, in_h, in_w, p_d, p_h, p_w);
          // if tmp_out_w > p_out_w  padding 1 more column 0s
          inputPaddingRightSide(prePaddingComp, inComp, tmp_out_d, tmp_out_h, p_out_w, 0, 0, tmp_out_w - p_out_w);
          for (int d = 0 ; d < tmp_out_d  ; d ++ ){
              for (int i = 0 ; i < tmp_out_h  ; i ++ ){
                  for (int j = 0 ; j < tmp_out_w ; j ++ ){      
                      std::cout << inComp[d * tmp_out_h * tmp_out_w +  i * tmp_out_w + j].real() << "," ; 
                  }
                   std::cout  << std::endl;        
                }
             std::cout  << std::endl;
          }
          std::cout  << std::endl;

          transReal2Comp(input_W + kernel_offset,keComp,tmp_out_d,tmp_out_h,tmp_out_w, k_d, k_h, k_w);

         // std::cout << "inComp" << std::endl;

          for (int d = 0 ; d < tmp_out_d  ; d ++ ){
              for (int i = 0 ; i < tmp_out_h  ; i ++ ){
                  for (int j = 0 ; j < tmp_out_w ; j ++ ){      
                    // std::cout << inComp[d * tmp_out_h * tmp_out_w +  i * tmp_out_w + j].real() << "," ; 
                  }
                  // std::cout  << std::endl;        
                }
            // std::cout  << std::endl;
          }

        //  std::cout << "keComp" << std::endl;
          for (int d = 0 ; d < tmp_out_d  ; d ++ ){
              for (int i = 0 ; i < tmp_out_h  ; i ++ ){
                  for (int j = 0 ; j < tmp_out_w ; j ++ ){        
                    // std::cout << keComp[d * tmp_out_h * tmp_out_w + i * tmp_out_w + j].real() << "," ; 
                  }
                  // std::cout  << std::endl;        
                }
              // std::cout  << std::endl; 
        }


          // fft a
          fft_cuda_3d(1,tmp_out_d,tmp_out_h,tmp_out_w,inComp,outinComp);
          // fft b
          fft_cuda_3d(1,tmp_out_d,tmp_out_h,tmp_out_w,keComp,outkeComp);
         
          oppositeImagPart(outkeComp,tmp_out_d * tmp_out_h * tmp_out_w);

          std::cout << "outinComp" << std::endl;

          for (int d = 0 ; d < tmp_out_d  ; d ++ ){
              for (int i = 0 ; i < tmp_out_h  ; i ++ ){
                  for (int j = 0 ; j < tmp_out_w ; j ++ ){      
                     std::cout << outinComp[d * tmp_out_h * tmp_out_w +  i * tmp_out_w + j].real() << "," ; 
                  }
                   std::cout  << std::endl;        
                }
             std::cout  << std::endl;
          }

          std::cout << "outkeComp" << std::endl;
          for (int d = 0 ; d < tmp_out_d  ; d ++ ){
              for (int i = 0 ; i < tmp_out_h  ; i ++ ){
                  for (int j = 0 ; j < tmp_out_w ; j ++ ){        
                     std::cout << outkeComp[d * tmp_out_h * tmp_out_w + i * tmp_out_w + j].real() << "," ; 
                  }
                   std::cout  << std::endl;        
                }
               std::cout  << std::endl; 
        }

          vectorCrossMultip(outinComp,outkeComp,inOutComp,tmp_out_d * tmp_out_h * tmp_out_w);
          // ifft c
          fft_cuda_3d(0,tmp_out_d,tmp_out_h,tmp_out_w, inOutComp, knkcOutOutComp + tmp_out_offset);



} 
        matrixAddition(knkcOutOutComp + tmp_out_kn_offset , k_c , tmp_out_d , tmp_out_h ,tmp_out_w);
        eliminateInputChannel(knkcOutOutComp + tmp_out_kn_offset, knOutOutComp + out_kn_offset, tmp_out_d, tmp_out_h, tmp_out_w);

        // std::cout << "knOutOutComp " << std::endl;
        for (int d = 0 ; d < tmp_out_d  ; d ++ ){
            for (int i = 0 ; i < tmp_out_h  ; i ++ ){
                for (int j = 0 ; j < tmp_out_w ; j ++ ){
                  // std::cout << knOutOutComp[ out_kn_offset + d * tmp_out_h * tmp_out_w + i * tmp_out_w + j].real() << "," ; 
                }
                // std::cout  << std::endl; 
            }
            // std::cout  << std::endl; 
        }
      }
    }
  //  std::cout << "all_knOutOutComp " << std::endl;
  //  std::cout << "params " << in_n << k_n << tmp_out_d << tmp_out_h << tmp_out_w << std::endl;

    for (int i = 0 ; i < in_n ; i ++){
      for (int j = 0 ; j < k_n ; j ++){
        for (int k = 0 ; k < tmp_out_d ; k ++){
          for (int l = 0 ; l < tmp_out_h ; l ++){
            for (int m = 0 ; m < tmp_out_w ; m ++){
              // std::cout << knOutOutComp[i * k_n * tmp_out_d * tmp_out_h * tmp_out_w + j * tmp_out_d * tmp_out_h * tmp_out_w + k * tmp_out_h * tmp_out_w + l * tmp_out_w + m].real() << "," ; 
            }
            // std::cout  << std::endl;
          }
          // std::cout  << std::endl; 
          }
        // std::cout  << std::endl; 
      }
      // std::cout  << std::endl; 
    }


    // tmp_out_w = p_out_w + 1;
    // 如果裁边tmp_out_w - p_out_w == 1 需要先去掉width轴最右侧的一部分
    Comp* knOutOutComp_trim_width = (Comp*) malloc( in_n * k_n * tmp_out_d * tmp_out_h * p_out_w * sizeof(Comp));
    matrixTrimRightWidth(knOutOutComp,knOutOutComp_trim_width, in_n , k_n , tmp_out_d , tmp_out_h ,tmp_out_w , p_out_w);
    // std::cout << "matrixTrimRightWidth" << std::endl;

    for (int i = 0 ; i < in_n ; i ++){
      for (int j = 0 ; j < k_n ; j ++){
        for (int k = 0 ; k < tmp_out_d ; k ++){
          for (int l = 0 ; l < tmp_out_h ; l ++){
            for (int m = 0 ; m < p_out_w ; m ++){
              // std::cout << knOutOutComp_trim_width[i * k_n * tmp_out_d * tmp_out_h * p_out_w + j * tmp_out_d * tmp_out_h * p_out_w + k * tmp_out_h * p_out_w + l * p_out_w + m].real() << "," ; 
            }
            // std::cout  << std::endl;
          }
          // std::cout  << std::endl; 
          }
        // std::cout  << std::endl; 
      }
      // std::cout  << std::endl; 
    }
    

    // 对输出结果进行裁边处理
    Comp* outOutTrimComp = (Comp*) malloc(out_n * out_c * out_d * out_h * out_w * sizeof(Comp));
    // Comp* outOutTrimComp = (Comp*) malloc(out_n * out_c * out_d * tmp_out_h * tmp_out_w * sizeof(Comp));
    matrixTrimPadding(knOutOutComp_trim_width, outOutTrimComp, out_n, out_c, tmp_out_d ,tmp_out_h, p_out_w, out_d ,out_h, out_w);


    // std::cout << "matrixTrimPadding " << std::endl;


    // for (int i = 0 ; i < out_n ; i ++){
    //   for (int j = 0 ; j < out_c ; j ++){
    //     for (int m = 0 ; m < out_d ; m ++){
    //         for (int k = 0 ; k < tmp_out_h ; k ++){
    //           for (int l = 0 ; l < tmp_out_w ; l ++){
    //             std::cout << outOutTrimComp[i * out_c * out_d * tmp_out_h * tmp_out_w + j * out_d * tmp_out_h * tmp_out_w + m * tmp_out_h * tmp_out_w + k * tmp_out_w + l].real() << "," ; 
    //           }
    //           std::cout  << std::endl; 
    //         }
    //     std::cout  << std::endl; 
    //   }
    //   std::cout  << std::endl; 
    // }
    //   std::cout  << std::endl; 
    // }



    for (int i = 0 ; i < out_n ; i ++){
      for (int j = 0 ; j < out_c ; j ++){
        for (int m = 0 ; m < out_d ; m ++){
            for (int k = 0 ; k < out_h ; k ++){
              for (int l = 0 ; l < out_w ; l ++){
                std::cout << outOutTrimComp[i * out_c * out_d * out_h * out_w + j * out_d * out_h * out_w + m * out_h * out_w + k * out_w + l].real() << "," ; 
              }
               std::cout  << std::endl; 
            }
         std::cout  << std::endl; 
      }
       std::cout  << std::endl; 
    }
     std::cout  << std::endl; 
    }



    // transform complex to real
    transComp2Real(outOutTrimComp, output_Y, out_n * out_d*  out_c * out_h * out_w);
}

void convfft_3d_v2(Real* input_X, int in_n, int in_c, int in_d, int in_h, int in_w,
  Real* input_W, int k_n, int k_c, int k_d, int k_h, int k_w,
  int p_d, int p_h, int p_w,
  int s_h, int s_w,
  int d_h, int d_w,
  Real* output_Y,int out_n, int out_c, int out_d, int out_h, int out_w)
{
    /**
      限制输入矩阵的通道数和kernel的输入通道数一致
    */
    if (k_c != in_c){
      throw string("convfft_3d cannot perform convolution because the input size does not match the kernel size!");
    }
   
    // tmp size of output
    size_t tmp_out_h = in_h + p_h * 2;
    size_t tmp_out_d = in_d + p_d * 2;
    size_t tmp_out_w = in_w + p_w * 2;


    tmp_out_h = pow(2,((int)log2(tmp_out_h - 1) + 1)); 
    tmp_out_w = pow(2,((int)log2(tmp_out_w - 1) + 1)); 
    tmp_out_d = pow(2,((int)log2(tmp_out_d - 1) + 1)); 

    
    Comp* knkcOutOutComp = (Comp*) malloc( in_n * k_n * k_c * tmp_out_d * tmp_out_h * tmp_out_w * sizeof(Comp));
    Comp* knOutOutComp = (Comp*) malloc( in_n * k_n * tmp_out_d * tmp_out_h * tmp_out_w * sizeof(Comp));
    Comp* prePaddingComp = (Comp*) malloc(in_n * in_c * (in_d + p_d * 2) * (in_h + p_h * 2) * (in_w + p_w * 2) * sizeof(Comp));

    Comp* inComp = (Comp*) malloc(in_n * in_c * tmp_out_d * tmp_out_h * tmp_out_w * sizeof(Comp));
    Comp* keComp = (Comp*) malloc(k_n * k_c * tmp_out_d * tmp_out_h * tmp_out_w * sizeof(Comp));

    Comp* outinComp = (Comp*) malloc(in_n * in_c * tmp_out_d * tmp_out_h * tmp_out_w * sizeof(Comp));
    Comp* outkeComp = (Comp*) malloc(k_n * k_c * tmp_out_d * tmp_out_h * tmp_out_w * sizeof(Comp));
    Comp* inOutComp = (Comp*) malloc(in_n * in_c * k_n * tmp_out_d * tmp_out_h * tmp_out_w * sizeof(Comp));

    
    memset(prePaddingComp,0,in_n * in_c * (in_d + p_d * 2) * (in_h + p_h * 2) * (in_w + p_w * 2) * sizeof(Comp));
    inputPadding(input_X, prePaddingComp, in_n, in_c, in_d, in_h, in_w , p_d, p_h, p_w);
    transReal2CompWithMemset(prePaddingComp, inComp, in_n, in_c, tmp_out_d, tmp_out_h, tmp_out_w, in_d + p_d * 2, in_h + p_h * 2, in_w + p_w * 2);
    transReal2CompWithMemset(input_W, keComp, k_n, k_c, tmp_out_d, tmp_out_h, tmp_out_w, k_d, k_h, k_w);
    
    
    cfft_plan plan_h = make_cfft_plan (tmp_out_h);
    
    cfft_plan plan_w = make_cfft_plan (tmp_out_w);
    
    cfft_plan plan_d = make_cfft_plan (tmp_out_d);
    
    // fft a
    for (int i = 0 ; i < in_n ; i++){
        for (int j = 0 ; j < in_c ; j++){
	    int offset = i * in_c * tmp_out_d * tmp_out_h * tmp_out_w + j * tmp_out_d * tmp_out_h * tmp_out_w ; 
    	    fft_cuda_3d_v2(plan_h, plan_w, plan_d, 1, tmp_out_d, tmp_out_h, tmp_out_w, inComp + offset , outinComp + offset);
	}
    }
    // fft b
    for (int i = 0 ; i < k_n ; i++){
        for (int j = 0 ; j < k_c ; j++){
	    int offset = i * k_c * tmp_out_d * tmp_out_h * tmp_out_w + j * tmp_out_d * tmp_out_h * tmp_out_w ; 
            fft_cuda_3d_v2(plan_h, plan_w, plan_d, 1,tmp_out_d,tmp_out_h,tmp_out_w,keComp + offset, outkeComp + offset);
	}
    }
    // give the image part an oppsite symbol
    oppositeImagPart(outkeComp, k_n * k_c * tmp_out_d * tmp_out_h * tmp_out_w);
   // vectorCrossMultip 
    for (int nn = 0 ; nn < in_n ; nn++){
        for (int k_index = 0 ; k_index < k_n ; k_index++){ 
            for (int c_index = 0 ; c_index < k_c ; c_index++){
                // the offset of each feature map
                int input_offset = nn * in_c * tmp_out_d * tmp_out_h * tmp_out_w + c_index  * tmp_out_d * tmp_out_h * tmp_out_w;
                int filter_offset = k_index * k_c * tmp_out_d * tmp_out_h * tmp_out_w + c_index * tmp_out_d * tmp_out_h * tmp_out_w;
                vectorCrossMultip(outinComp + input_offset, outkeComp + filter_offset, inOutComp + nn * k_n * k_c * tmp_out_d * tmp_out_h * tmp_out_w + k_index * k_c * tmp_out_d * tmp_out_h * tmp_out_w + c_index * tmp_out_d * tmp_out_h * tmp_out_w, tmp_out_d * tmp_out_h * tmp_out_w);
	    }
	}
    }
    // fft 2d inverse
    for (int nn = 0 ; nn < in_n ; nn++){
        for (int k_index = 0 ; k_index < k_n ; k_index++){ 
            for (int c_index = 0 ; c_index < k_c ; c_index++){
                // the offset of each feature map
                fft_cuda_3d_v2(plan_h, plan_w, plan_d, 0, tmp_out_d, tmp_out_h, tmp_out_w, 
				inOutComp+ nn * k_n * k_c * tmp_out_d * tmp_out_h * tmp_out_w + k_index * k_c * tmp_out_d * tmp_out_h * tmp_out_w + c_index * tmp_out_d  * tmp_out_h * tmp_out_w, 
				knkcOutOutComp + nn * k_n * k_c * tmp_out_d  * tmp_out_h * tmp_out_w + k_index * k_c * tmp_out_d * tmp_out_h * tmp_out_w + c_index* tmp_out_d * tmp_out_h * tmp_out_w );
	    }
	}
    }
   // eliminate input channel 
   for (int nn = 0 ; nn < in_n ; nn++){
        int n_kn_kc_out_offset = nn * k_n * k_c * tmp_out_d * tmp_out_h * tmp_out_w;
        // the output channel of kernel filter
        for (int k_index = 0 ; k_index < k_n ; k_index++){ 
            int tmp_out_kn_offset = n_kn_kc_out_offset + k_index * k_c * tmp_out_d * tmp_out_h * tmp_out_w;
            int out_kn_offset = nn * k_n * tmp_out_d * tmp_out_h * tmp_out_w + k_index * tmp_out_d * tmp_out_h * tmp_out_w;
            matrixAddition(knkcOutOutComp + tmp_out_kn_offset , k_c, tmp_out_d, tmp_out_h, tmp_out_w);
            eliminateInputChannel(knkcOutOutComp + tmp_out_kn_offset, knOutOutComp + out_kn_offset, tmp_out_d, tmp_out_h, tmp_out_w);
        }
    }
    
    // 对输出结果进行裁边处理
    Comp* outOutTrimComp = (Comp*) malloc(out_n * out_c * out_d * out_h * out_w * sizeof(Comp));
    matrixTrimPadding(knOutOutComp, outOutTrimComp, out_n, out_c, tmp_out_d, tmp_out_h, tmp_out_w, out_d, out_h, out_w);
    // transform complex to real
    transComp2Real(outOutTrimComp, output_Y, out_n * out_c * out_d * out_h * out_w);
}
// this is the version which just to fit the backdata_3d 
// try to change the position of zero padding until get the final result fit to the cudnn 
void convfft_3d_4backdata(Real* input_X, int in_n, int in_c, int in_d, int in_h, int in_w,
  Real* input_W, int k_n, int k_c, int k_d, int k_h, int k_w,
  int p_d, int p_h, int p_w,
  int s_h, int s_w,
  int d_h, int d_w,
  Real* output_Y,int out_n, int out_c, int out_d, int out_h, int out_w)
{
    /**
      限制输入矩阵的通道数和kernel的输入通道数一致
    */
    if (k_c != in_c){
      cout << "convfft_3d cannot perform convolution because the input size does not match the kernel size!" << endl;
      throw string("convfft_3d cannot perform convolution because the input size does not match the kernel size!");
    }
   
    // tmp size of output
    // size_t tmp_out_h = in_h + p_h * 2;
    // 
    // // size_t tmp_out_d = in_d + k_d - 1;
    // size_t tmp_out_d = in_d + p_d * 2;
    // size_t tmp_out_w = in_w + p_w * 2;
    size_t tmp_out_h = in_h + k_h - 1;
    size_t tmp_out_w = in_w + k_w - 1;
    size_t tmp_out_d = in_d + k_d - 1;

    std::cout << "tmp_out_h: "<< tmp_out_h << std::endl;  
    std::cout << "tmp_out_w: "<< tmp_out_w << std::endl; 
    std::cout << "tmp_out_d: "<< tmp_out_d << std::endl; 

    tmp_out_h = pow(2,((int)log2(tmp_out_h - 1) + 1)); 
    tmp_out_w = pow(2,((int)log2(tmp_out_w - 1) + 1)); 
    tmp_out_d = pow(2,((int)log2(tmp_out_d - 1) + 1)); 

    std::cout << "tmp_out_h: "<< tmp_out_h << std::endl;  
    std::cout << "tmp_out_w: "<< tmp_out_w << std::endl; 
    std::cout << "tmp_out_d: "<< tmp_out_d << std::endl; 
    
    Comp* knkcOutOutComp = (Comp*) malloc( in_n * k_n * k_c * tmp_out_d * tmp_out_h * tmp_out_w * sizeof(Comp));
    Comp* knOutOutComp = (Comp*) malloc( in_n * k_n * tmp_out_d * tmp_out_h * tmp_out_w * sizeof(Comp));
    Comp* prePaddingComp = (Comp*) malloc(in_n * in_c * (in_d + p_d * 2) * (in_h + p_h * 2) * (in_w + p_w * 2) * sizeof(Comp));

    Comp* inComp = (Comp*) malloc(in_n * in_c * tmp_out_d * tmp_out_h * tmp_out_w * sizeof(Comp));
    Comp* keComp = (Comp*) malloc(k_n * k_c * tmp_out_d * tmp_out_h * tmp_out_w * sizeof(Comp));

    Comp* outinComp = (Comp*) malloc(in_n * in_c * tmp_out_d * tmp_out_h * tmp_out_w * sizeof(Comp));
    Comp* outkeComp = (Comp*) malloc(k_n * k_c * tmp_out_d * tmp_out_h * tmp_out_w * sizeof(Comp));
    Comp* inOutComp = (Comp*) malloc(in_n * in_c * k_n * tmp_out_d * tmp_out_h * tmp_out_w * sizeof(Comp));

    
    memset(prePaddingComp,0,in_n * in_c * (in_d + p_d * 2) * (in_h + p_h * 2) * (in_w + p_w * 2) * sizeof(Comp));
    inputPadding(input_X, prePaddingComp, in_n, in_c, in_d, in_h, in_w , p_d, p_h, p_w);

    transReal2CompWithMemset(prePaddingComp, inComp, in_n, in_c, tmp_out_d, tmp_out_h, tmp_out_w, in_d + p_d * 2, in_h + p_h * 2, in_w + p_w * 2);
    transReal2CompWithMemset(input_W, keComp, k_n, k_c, tmp_out_d, tmp_out_h, tmp_out_w, k_d, k_h, k_w);
    
   cout << "inComp" << endl;
   for (int nn = 0 ; nn < in_n ; nn++){
            for (int c_index = 0 ; c_index < in_c ; c_index++){
                // the offset of each feature map
               for (int k = 0 ; k < tmp_out_d ; k++){
                   for (int i = 0 ; i < tmp_out_h ; i++){
		      for (int j = 0 ; j < tmp_out_w ; j++){
		          cout << (inComp + nn * in_c * tmp_out_d * tmp_out_h * tmp_out_w + c_index * tmp_out_d * tmp_out_h * tmp_out_w )[ k * tmp_out_h * tmp_out_w +  i * tmp_out_w + j].real()<<",";
		      }
		      cout << endl;
		   }
		      cout << endl;
	       }
	      cout << endl;
	  }
        cout << endl;
    }
   cout << endl;
    
    cout << "make_cfft_plan plan_h" << endl;
    cfft_plan plan_h = make_cfft_plan (tmp_out_h);
    
    cout << "make_cfft_plan plan_w" << endl;
    cfft_plan plan_w = make_cfft_plan (tmp_out_w);
    
    cout << "make_cfft_plan plan_d" << endl;
    cfft_plan plan_d = make_cfft_plan (tmp_out_d);
    
    cout << "fft_cuda_3d_v2 for input" << endl;
    // fft a
    for (int i = 0 ; i < in_n ; i++){
        for (int j = 0 ; j < in_c ; j++){
	    int offset = i * in_c * tmp_out_d * tmp_out_h * tmp_out_w + j * tmp_out_d * tmp_out_h * tmp_out_w ; 
    	    fft_cuda_3d_v2(plan_h, plan_w, plan_d, 1, tmp_out_d, tmp_out_h, tmp_out_w, inComp + offset , outinComp + offset);
	}
    }
   cout << "outinComp" << endl;
   for (int nn = 0 ; nn < in_n ; nn++){
            for (int c_index = 0 ; c_index < in_c ; c_index++){
                // the offset of each feature map
               for (int k = 0 ; k < tmp_out_d ; k++){
                   for (int i = 0 ; i < tmp_out_h ; i++){
		      for (int j = 0 ; j < tmp_out_w ; j++){
		          cout << (outinComp + nn * in_c * tmp_out_d * tmp_out_h * tmp_out_w + c_index * tmp_out_d * tmp_out_h * tmp_out_w )[ k * tmp_out_h * tmp_out_w +  i * tmp_out_w + j].real()<<",";
		      }
		      cout << endl;
		   }
		      cout << endl;
	       }
	      cout << endl;
	  }
        cout << endl;
    }
   cout << endl;

    cout << "fft_cuda_3d_v2 for weight" << endl;
    // fft b
    for (int i = 0 ; i < k_n ; i++){
        for (int j = 0 ; j < k_c ; j++){
	    int offset = i * k_c * tmp_out_d * tmp_out_h * tmp_out_w + j * tmp_out_d * tmp_out_h * tmp_out_w ; 
            fft_cuda_3d_v2(plan_h, plan_w, plan_d, 1,tmp_out_d,tmp_out_h,tmp_out_w,keComp + offset, outkeComp + offset);
	}
    }
    cout << "oppositeImagPart for weight" << endl;
    // give the image part an oppsite symbol
    oppositeImagPart(outkeComp, k_n * k_c * tmp_out_d * tmp_out_h * tmp_out_w);
    
    cout << "vectorCrossMultip for weight" << endl;
   // vectorCrossMultip 
    for (int nn = 0 ; nn < in_n ; nn++){
        for (int k_index = 0 ; k_index < k_n ; k_index++){ 
            for (int c_index = 0 ; c_index < k_c ; c_index++){
                // the offset of each feature map
                int input_offset = nn * in_c * tmp_out_d * tmp_out_h * tmp_out_w + c_index  * tmp_out_d * tmp_out_h * tmp_out_w;
                int filter_offset = k_index * k_c * tmp_out_d * tmp_out_h * tmp_out_w + c_index * tmp_out_d * tmp_out_h * tmp_out_w;
                vectorCrossMultip(outinComp + input_offset, outkeComp + filter_offset, inOutComp + nn * k_n * k_c * tmp_out_d * tmp_out_h * tmp_out_w + k_index * k_c * tmp_out_d * tmp_out_h * tmp_out_w + c_index * tmp_out_d * tmp_out_h * tmp_out_w, tmp_out_d * tmp_out_h * tmp_out_w);
	    }
	}
    }
   cout << "inOutComp real part is " << endl;
   for (int nn = 0 ; nn < in_n ; nn++){
        for (int k_index = 0 ; k_index < k_n ; k_index++){ 
            for (int c_index = 0 ; c_index < k_c ; c_index++){
                // the offset of each feature map
               for (int k = 0 ; k < tmp_out_d ; k++){
               for (int i = 0 ; i < tmp_out_h ; i++){
		   for (int j = 0 ; j < tmp_out_w ; j++){
		       cout << (inOutComp+ nn * k_n * k_c * tmp_out_d * tmp_out_h * tmp_out_w + k_index * k_c * tmp_out_d * tmp_out_h * tmp_out_w + c_index* tmp_out_d* tmp_out_h * tmp_out_w )[k * tmp_out_h * tmp_out_w + i * tmp_out_w + j].real()<<",";
		   }
		   cout << endl;
		}

		   cout << endl;
	    }
		   cout << endl;
	    }
		   cout << endl;
	}
		   cout << endl;
    }
		   cout << endl;
    cout << "fft_cuda_3d_v2 inverse" << endl;
    // fft 2d inverse
    for (int nn = 0 ; nn < in_n ; nn++){
        for (int k_index = 0 ; k_index < k_n ; k_index++){ 
            for (int c_index = 0 ; c_index < k_c ; c_index++){
                // the offset of each feature map
                fft_cuda_3d_v2(plan_h, plan_w, plan_d, 0, tmp_out_d, tmp_out_h, tmp_out_w, 
				inOutComp+ nn * k_n * k_c * tmp_out_d * tmp_out_h * tmp_out_w + k_index * k_c * tmp_out_d * tmp_out_h * tmp_out_w + c_index * tmp_out_d  * tmp_out_h * tmp_out_w, 
				knkcOutOutComp + nn * k_n * k_c * tmp_out_d  * tmp_out_h * tmp_out_w + k_index * k_c * tmp_out_d * tmp_out_h * tmp_out_w + c_index* tmp_out_d * tmp_out_h * tmp_out_w );
	    }
	}
    }
   cout << "knkcOutOutComp" << endl;
   for (int nn = 0 ; nn < in_n ; nn++){
        for (int k_index = 0 ; k_index < k_n ; k_index++){ 
            for (int c_index = 0 ; c_index < k_c ; c_index++){
                // the offset of each feature map
               for (int k = 0 ; k < tmp_out_d; k++){
               for (int i = 0 ; i < tmp_out_h ; i++){
		   for (int j = 0 ; j < tmp_out_w ; j++){
		       cout << (knkcOutOutComp+ nn * k_n * k_c * tmp_out_d  * tmp_out_h * tmp_out_w + k_index * k_c * tmp_out_d * tmp_out_h * tmp_out_w + c_index * tmp_out_d * tmp_out_h * tmp_out_w )[k * tmp_out_h * tmp_out_w + i * tmp_out_w + j].real()<<",";
		   }
		   cout << endl;
		}
		   cout << endl;
	       }
		   cout << endl;
	    }
		   cout << endl;
	}
		   cout << endl;
    }
		   cout << endl;

    cout << "eliminateInputChannel" << endl;
   // eliminate input channel 
   for (int nn = 0 ; nn < in_n ; nn++){
        int n_kn_kc_out_offset = nn * k_n * k_c * tmp_out_d * tmp_out_h * tmp_out_w;
        // the output channel of kernel filter
        for (int k_index = 0 ; k_index < k_n ; k_index++){ 
            int tmp_out_kn_offset = n_kn_kc_out_offset + k_index * k_c * tmp_out_d * tmp_out_h * tmp_out_w;
            int out_kn_offset = nn * k_n * tmp_out_d * tmp_out_h * tmp_out_w + k_index * tmp_out_d * tmp_out_h * tmp_out_w;
            matrixAddition(knkcOutOutComp + tmp_out_kn_offset , k_c, tmp_out_d, tmp_out_h, tmp_out_w);
            eliminateInputChannel(knkcOutOutComp + tmp_out_kn_offset, knOutOutComp + out_kn_offset, tmp_out_d, tmp_out_h, tmp_out_w);
        }
    }
    
    cout << "matrixTrimPadding" << endl;
    // 对输出结果进行裁边处理
    Comp* outOutTrimComp = (Comp*) malloc(out_n * out_c * out_d * out_h * out_w * sizeof(Comp));
    matrixTrimPadding(knOutOutComp, outOutTrimComp, out_n, out_c, tmp_out_d, tmp_out_h, tmp_out_w, out_d, out_h, out_w);
    for (int i = 0 ; i < out_n ; i ++){
      for (int j = 0 ; j < out_c ; j ++){
        for (int m = 0 ; m < out_d ; m ++){
            for (int k = 0 ; k < out_h ; k ++){
              for (int l = 0 ; l < out_w ; l ++){
                std::cout << outOutTrimComp[i * out_c * out_d * out_h * out_w + j * out_d * out_h * out_w + m * out_h * out_w + k * out_w + l].real() << "," ; 
              }
               std::cout  << std::endl; 
            }
         std::cout  << std::endl; 
      }
       std::cout  << std::endl; 
    }
     std::cout  << std::endl; 
    }
    // transform complex to real
    transComp2Real(outOutTrimComp, output_Y, out_n * out_c * out_d * out_h * out_w);
}

#endif
